# ===============================================================================================================================
# ===============================================================================================================================
# ACS Analysis
# This code is for replicating: Appendix Table A3
# 
#   1. Import and Prep ACS data
#   2. Basic analysis
#   2. Create Table
# ===============================================================================================================================
# ===============================================================================================================================


#Clear out the namespace
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects. (simpler is rm(list=ls()))
gc() #free up memory and report the memory usage.
.pardefault <- par()


#Turn off scientific notation (http://datacornering.com/remove-scientific-notation-in-r/)
# options(scipen = 100)
options(scipen = 999)

# Set the path to the root folder ending in RFSDelivery 
cd <- "D:/..."
cd <- "D:/jodo_emga/RFSDelivery/"
# Initialize all other paths and import all functions and packages from the AllSettings file.
source(paste0(cd,'/Code/R/AllSettings.R',collapse=''))


# ===============================================================================================================================
# ===============================================================================================================================
# Import the ACS dataset
# ===============================================================================================================================
# ===============================================================================================================================
acs <- read_parquet(cdd['p_d_acs_pa'] %c% '/acs_zip_analysis.parquet',sep='')

acs <- acs %>% mutate(
  pop_l = log(population),
  inc_hhmean_l = log(inc_hhmean),
  mhi_quin1_l = log(mhi_quin1)
)

# Variable labels
labels_acs = list(zip='Zip Code',year='Year',population='Population',pop_l='Log Population',male_p='Male',white_p='White',
                  educ_sc_p='Some College',educ_bach_p='Bachelors',pvf100_p='Poor <100FPL',pvf200_p='Poor <200FPL',
                  ownocc_p='Homeowner',gini='Gini',insured_p='Insured',snap_p='SNAP',aid_p='Aid',pubai_p='Public Asst.',
                  inc_hhmean='HH Inc Mean',inc_hhmean_l='Log HH Inc Mean',mhi_quin1='HH Inc 1st Quin',mhi_quin1_l='Log HH Inc 1st Quin',
                  comm_public_p='Comm Transit',comm_car_p='Comm Car',comm_timemean='Comm Time',
                  empt_nw_p='Not Working', empt_nft_p='Working Part Time', empt_ft_p='Working Full Time')
labels_coh = list(treated='Treated',etime='Event Time YQ',etimey='Event Time Y',distance_t2cf='Dist of Target PC to CF PC',date='Date',
                  ctype='Counterfactual type',cohort='PC Cohort')
labels_indep = list(distance='Dist',ic='IC')
labels_other = list(pop2sqkm_wq10='Pop Dens Dec')
labels_acs <- c(labels_acs, labels_coh, labels_indep, labels_other)

# Apply variable labels
labels_acs_sub <- labels_acs[names(labels_acs) %in% names(acs)]
var_label(acs) <- labels_acs_sub



#Set the labels used for all tables
setFixest_dict(unlist(labels_acs_sub))
my.style.df = style.df(fixef.title = "", fixef.suffix = " FE", stats.title = "", slopes.title = "")
setFixest_etable(digits = "r4", digits.stats = "r4", style.df = my.style.df)


#Turn graphics off to allow saving the plots
graphics.off()



# ===============================================================================================================================
# Estimate and Output Appendix Table A3
# ===============================================================================================================================
spec_y = list(
  'table_vars'=c('pop_l','pvf100_p','empt_nw_p', 'empt_nft_p', 'empt_ft_p', 'mhi_quin1_l', 'pubai_p')
)
spec_fe = list(
  'cz_ctpd10_sd'='cohort^zip + cohort^etimey^pop2sqkm_wq10 + state^date',
  'cz_ct_sd'='cohort^zip + cohort^etimey + state^date'
)
spec_treat <- list('Tty' = 'i(etimey, treated, ref = -1)')
var_i <- 'Tty'
spec_trl <- 'l'
var_tr <- 'd50_ic50'

form_y <- 'c(' %c% paste0(unlist(spec_y[['table_vars']]),collapse=', ') %c% ') '
form_ife <- spec_treat[var_i] %c% ' | ' %c% spec_fe['cz_ct_sd']
formula_est_base <- as.formula(form_y %c% ' ~ ' %c% form_ife)
est_base <- feols(formula_est_base, cluster=~cohort, acs)

form_y <- 'c(' %c% paste0(unlist(spec_y['table_vars']),collapse=', ') %c% ') '
form_ife <- spec_treat[var_i] %c% ' | ' %c% spec_fe['cz_ctpd10_sd']
formula_est_aug <- as.formula(form_y %c% ' ~ ' %c% form_ife)
est_aug <- feols(formula_est_aug, cluster=~cohort, acs)


#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Appendix Table A3: Demographic trends in treated versus control zip codes, annual stacked DiD estimates
summary_table <- etable(c(as.list(est_base),as.list(est_aug)))
etable(c(as.list(est_base),as.list(est_aug)), style.tex = style.tex("qje"), 
       signif.code=c("***"=0.01,"**"=0.05,"*"=0.10), fitstat = ~ r2 + n, tex = TRUE, 
       file = cdd['p_rt'] %c% '/General/Appendix Table A3.tex', replace = TRUE)









